{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c0912d9e-45db-4cc9-9880-7f2c44cfd369",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Case4 import main\n",
    "from fem import Poisson2D\n",
    "import calfem.geometry as cfg\n",
    "import calfem.mesh as cfm\n",
    "import calfem.vis as cfv\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2e6b12fa-c62b-41b3-ac30-49c921fa6f14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# def create_geometry(el_on_curve=20):\n",
    "#     # externel boundary ID\n",
    "#     dir_id_1 = 100\n",
    "#     dir_id_2 = 200\n",
    "#     dir_id_3 = 300\n",
    "#     neu_id_1 = 400\n",
    "#     # hole boundaryd ID\n",
    "#     hole_neu_id_1 = 500\n",
    "#     hole_neu_id_2 = 600\n",
    "#     hole_neu_id_3 = 700\n",
    "#     hole_neu_id_4 = 800\n",
    "\n",
    "#     # Geometry definition\n",
    "#     g = cfg.Geometry()\n",
    "#     g.point([0.0, 0.0], ID=0)\n",
    "#     g.point([1.0, 0.0], ID=1)\n",
    "#     g.point([1.0, 1.0], ID=2)\n",
    "#     g.point([0.0, 1.0], ID=3)\n",
    "#     g.point([0.4, 0.4], ID=4)\n",
    "#     g.point([0.6, 0.4], ID=5)\n",
    "#     g.point([0.6, 0.6], ID=6)\n",
    "#     g.point([0.4, 0.6], ID=7)\n",
    "#     g.spline([0, 1], el_on_curve=el_on_curve, marker=100, ID=0)\n",
    "#     g.spline([1, 2], el_on_curve=el_on_curve, marker=101, ID=1) \n",
    "#     g.spline([2, 3], el_on_curve=el_on_curve, marker=200, ID=2) \n",
    "#     g.spline([3, 0], el_on_curve=el_on_curve, marker=300, ID=3)\n",
    "#     g.spline([4, 5], el_on_curve=el_on_curve/4, marker=201, ID=4)\n",
    "#     g.spline([5, 6], el_on_curve=el_on_curve/4, marker=301, ID=5)\n",
    "#     g.spline([6, 7], el_on_curve=el_on_curve/4, marker=401, ID=6)\n",
    "#     g.spline([7, 4], el_on_curve=el_on_curve/4, marker=501, ID=7)\n",
    "#     g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "#     # cfv.drawGeometry(g)\n",
    "#     # cfv.showAndWait()\n",
    "\n",
    "\n",
    "#     bcs = {100:[0],\n",
    "#            200:[2],\n",
    "#            300:[3],\n",
    "#            101:[1],\n",
    "#            201:[4],\n",
    "#            301:[5],\n",
    "#            401:[6],\n",
    "#            501:[7],\n",
    "#            }\n",
    "\n",
    "#     for marker in bcs:\n",
    "#         for i in bcs[marker]:\n",
    "#             g.curve_marker(ID=i, marker=marker)\n",
    "#     # cfv.drawGeometry(g)\n",
    "#     # cfv.showAndWait()\n",
    "#     return g\n",
    "\n",
    "\n",
    "# def dir100(x,y):\n",
    "#     return np.zeros_like(x)\n",
    "\n",
    "# def dir200(x,y):\n",
    "#     return np.sin(x)/np.sin(1)\n",
    "\n",
    "# def dir300(x,y):\n",
    "#     return np.zeros_like(x)\n",
    "\n",
    "# def neu101(x,y):\n",
    "#     return (1/np.tan(1))*np.sinh(y)/np.sinh(1)\n",
    "\n",
    "# def neu201(x,y):\n",
    "#     return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "# def neu301(x,y):\n",
    "#     return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "# def neu401(x,y):\n",
    "#     return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "# def neu501(x,y):\n",
    "#     return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "\n",
    "# geometry = create_geometry(el_on_curve=20)\n",
    "# Dirichlet = [100,200,300]\n",
    "# Neumann = [101,201,301,401,501]\n",
    "# hole = [[201,301,401,501]]\n",
    "\n",
    "# BCfunc = {100:dir100, 200:dir200, 300:dir300,\n",
    "#           101:neu101, 201:neu201, 301:neu301, 401:neu401,501:neu501}\n",
    "\n",
    "\n",
    "# model = Poisson2D(geometry, Dirichlet, Neumann, hole, BCfunc)\n",
    "# model.mesh(showMesh=False)\n",
    "# model.efem()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "48c087bc-b40a-4b9d-9997-0cec2e87879b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# main()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c54063d-1519-4033-a57d-680bfe2cf388",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n",
      "size: 661 x 661\n",
      "rank: 660\n",
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x000001E9E6832520>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def create_geometry(el_on_curve=5):\n",
    "    # Geometry definition\n",
    "    g = cfg.Geometry()\n",
    "\n",
    "    points = [[0,1],[0.4,1], [0.6,1], [1,1],\n",
    "              [0,0.6], [0.4,0.6], [0.6,0.6], [1,0.6],\n",
    "              [0,0.4],[0.4,0.4],[0.6,0.4],[1,0.4],\n",
    "              [0,0], [0.4,0], [0.6,0], [1,0]]\n",
    "\n",
    "    splines = [[0,1],[1,2],[2,3], \n",
    "               [0,4],[1,5],[2,6],[3,7],\n",
    "               [4,5],[5,6],[6,7],\n",
    "               [4,8],[5,9],[6,10],[7,11],\n",
    "               [8,9],[9,10],[10,11],\n",
    "               [8,12],[9,13],[10,14],[11,15],\n",
    "               [12,13],[13,14],[14,15]\n",
    "              ]\n",
    "\n",
    "    longs = [0,2,3,4,5,6,7,9,14,16,17,18,19,20,21,23]\n",
    "\n",
    "    surfaces = [[0,4,7,3],[1,5,8,4],[2,6,9,5],[7,11,14,10],[9,13,16,12],[14,18,21,17],[15,19,22,18],\n",
    "             [16,20,23,19]]\n",
    "\n",
    "    for x, y in points:\n",
    "        g.point([x, y])\n",
    "\n",
    "    for i, s in enumerate(splines):\n",
    "        if i in longs:\n",
    "            g.spline(s, el_on_curve=2*el_on_curve)\n",
    "        else:\n",
    "            g.spline(s, el_on_curve=el_on_curve)\n",
    "\n",
    "    for f in surfaces:\n",
    "        g.struct_surf(f)\n",
    "\n",
    "\n",
    "    bcs = {100:[21,22,23],\n",
    "           200:[0,1,2],\n",
    "           300:[3,10,17],\n",
    "           101:[6,13,20],\n",
    "           201:[15],\n",
    "           301:[12],\n",
    "           401:[8],\n",
    "           501:[11],\n",
    "           }\n",
    "\n",
    "    for marker in bcs:\n",
    "        for i in bcs[marker]:\n",
    "            g.curve_marker(ID=i, marker=marker)\n",
    "    cfv.drawGeometry(g)\n",
    "    cfv.showAndWait()\n",
    "    return g\n",
    "\n",
    "\n",
    "def dir100(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir200(x,y):\n",
    "    return np.sin(x)/np.sin(1)\n",
    "\n",
    "def dir300(x,y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def neu101(x,y):\n",
    "    return (1/np.tan(1))*np.sinh(y)/np.sinh(1)\n",
    "\n",
    "def neu201(x,y):\n",
    "    return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "def neu301(x,y):\n",
    "    return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "def neu401(x,y):\n",
    "    return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "def neu501(x,y):\n",
    "    return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "geometry = create_geometry(el_on_curve=5)\n",
    "Dirichlet = [100,200,300]\n",
    "Neumann = [101,201,301,401,501]\n",
    "hole = [[201,301,401,501]]\n",
    "\n",
    "BCfunc = {100:dir100, 200:dir200, 300:dir300,\n",
    "          101:neu101, 201:neu201, 301:neu301, 401:neu401,501:neu501}\n",
    "\n",
    "\n",
    "model = Poisson2D(geometry, Dirichlet, Neumann, hole, BCfunc)\n",
    "model.mesh(showMesh=False)\n",
    "model.efem()\n",
    "\n",
    "def exact_func(x,y):\n",
    "    return (np.cos(x)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "model.display_exact(exact_func)\n",
    "# model.fem()\n",
    "\n",
    "# model.display()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d42dc183-1422-4277-9e3a-76e0c2bf25aa",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EFEM",
   "language": "python",
   "name": "efem"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
